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Abstract — In  recent  years  there  has  been  an  increase  in  the 
number  of  inactive  and  debris  objects  in  space.  The  characteri¬ 
zation  of  the  uncertainty  in  the  knowledge  of  these  Space  Objects 
(SOs)  is  very  important  in  developing  an  understanding  of  the 
space  debris  fields  and  any  present  or  future  threat  they  may 
pose.  This  work  examines  classification  based  on  Multiple  Model 
Adaptive  Estimation  (MMAE)  to  extract  SO  characteristics  from 
observations  while  estimating  the  probability  the  observations 
belong  to  a  given  class  of  objects.  Recovering  these  characteristics 
and  trajectories  with  sufficient  accuracy  is  shown  in  this  paper, 
where  the  characteristics  are  inherent  in  unique  SO  models  used 
in  the  MMAE  filter  bank.  A  number  of  scenarios  are  shown  to 
highlight  the  effectiveness  of  the  proposed  classification  approach. 
The  performance  of  this  strategy  is  demonstrated  via  simulated 
scenarios. 

I.  Introduction 

Due  to  the  large  number  of  space  objects  (SOs)  and  the 
limited  number  of  sensors  available  to  track  them,  it  is  difficult 
to  maintain  persistent  surveillance,  and,  therefore,  there  is 
inherent  uncertainty  and  latency  in  the  knowledge  of  the 
SO  population.  Although  the  amount  of  light  collected  from 
these  objects  is  small,  information  can  still  be  extracted  from 
photometric  data  which  can  be  used  to  determine  shapes  and 
other  properties.  Light  curve  data  are  the  time- varying  sensor 
wavelength-dependent  apparent  magnitude  of  energy  (e.g.  pho¬ 
tons)  scattered  (reflected)  off  of  an  object  along  the  line-of- 
sight  to  an  observer.  Attitude  estimation  and  extract  of  other 
characteristic  using  light  curve  data  has  been  demonstrated  in 
Refs.  1-6. 

An  approach  is  presented  that  uses  the  probability  from 
a  Multiple  Model  Adaptive  Estimation  (MMAE)  process  to 
determine  the  probability  that  a  given  Space  Object  (SO)  falls 
in  a  given  class.  MMAE  is  a  recursive  algorithm  that  uses  a 
bank  of  estimators,  each  dependent  on  a  particular  hypothesis, 
to  determine  an  estimate  based  upon  an  unknown  physical 
process  under  consideration.  In  particular,  the  hypotheses  can 
correspond  to  different  mathematical  models  of  the  same 
physical  process  or  of  the  same  model  but  dependent  upon 
different  constants  or  model  parameters.  The  classification 
approach  used  in  the  work  in  outlined  in  Figure  1.  The  first 
determination  is  made  from  the  size  of  the  shape  models  in  the 
bank.  For  each  model  in  the  bank  an  aspect  ratio  is  calculated 


for  each  size  by  calculating  the  length  ratio  of  that  side  with 
respect  to  the  largest  side  and  if  a  given  model  has  an  aspect 
ratio  less  than  0.1  it  is  considered  to  be  a  fragment. 


Fig.  1:  Overview  of  Classification  Approach 


The  second  classification  determination  is  made  from  the 
control  states.  For  each  model  in  the  bank  that  is  not  a  frag¬ 
ment  or  a  rocket  body  additional  models  are  created  that  are 
copies  of  that  shape  model  but  have  different  control  profiles. 
The  control  profiles  include  uncontrolled,  Sun  pointing,  spin 
stabilized,  and  nadir  pointing.  Therefore,  for  models  that  are 
not  fragments  of  rocket  bodies  3  additional  modes  are  created 
with  control.  The  control  states  are  not  limited  to  ones  used 
in  this  work  and  large  list  of  control  state  may  be  possible  but 
for  this  study  these  three  are  sufficient. 

The  classification  is  determined  using  the  shape  model,  for 
example  determining  whether  an  object  is  intact  or  passive  and 
whether  it  is  a  rocket  body  or  a  payload.  The  final  classification 
further  separates  the  SO  into  the  type  of  control  state  and 
determining  whether  an  object  is  uncontrolled  or  Sun  pointing, 
etc.  This  method  uses  the  MMAE  probability  to  classify  the 
four  feature  classes  and  results  for  this  method  are  shown.  This 
paper  discusses  the  theory  involved  behind  the  algorithm  and 
results  from  a  variety  of  simulation  trials  are  shown. 

The  organization  of  this  paper  is  as  follows.  First,  the 


Ashikmin- Shirley  light  curve  model  is  shown,  and  the  kine¬ 
matic  and  dynamic  models  used  in  this  work  are  discussed. 
Next,  the  Unscented  Kalman  Filter  approach  used  in  this  work 
is  outlined.  Following  this,  the  MMAE  approach  used  in  this 
work  is  outlined.  Then  the  classification  approach  is  discussed. 
Finally,  results  are  shown  for  simulated  examples.  Discussions 
and  conclusions  are  provided. 


II.  Ashikhmin-Shirley  Model 

Figure  2  shows  the  space  object  shape  model  and  reflection 
geometry.  In  addition  to  the  azimuth  and  elevation,  the  optical 
site  also  records  the  magnitude  of  the  brightness  of  the  Space 
Objects  (SOs).  The  brightness  of  an  object  in  space  can  be 
modeled  using  an  anisotropic  Phong  light  diffusion  model  or 
the  Ashikhmin-Shirley  model.7  This  model  is  based  on  the 
bidirectional  reflectance  distribution  function  (BRDF)  which 
models  light  distribution  scattered  from  the  surface  due  to 
the  incident  light.  The  BRDF  at  any  point  on  the  surface 
is  a  function  of  two  directions,  the  direction  from  which 
the  light  source  originates  and  the  direction  from  which  the 
scattered  light  leaves  the  observed  surface.  The  model  in  Ref.  7 
decomposes  the  BRDF  into  a  specular  component  and  a  diffuse 
component.  The  two  terms  sum  to  give  the  total  BRDF: 

Ptotal(^)  =  Pspec(^)  T~  Pdiff(^)  (1) 

where  i  denotes  the  ith  facet  of  the  SOs.  Each  facet  contributes 
independently  to  the  brightness  and  total  brightness  is  the 
sum  over  each  facet’s  contribution.  The  diffuse  component 
represents  light  that  is  scattered  equally  in  all  directions 
(Lambertian)  and  the  specular  component  represents  light  that 
is  concentrated  about  some  direction  (mirror- like).  Reference 
7  develops  a  model  for  continuous  arbitrary  surfaces  but 
simplifies  for  flat  surfaces.  This  simplified  model  is  employed 
in  this  work  as  shape  models  are  considered  to  consist  of 
a  finite  number  of  flat  facets.  Therefore  the  total  observed 
brightness  of  an  object  becomes  the  sum  of  the  contribution 
from  each  facet.  Under  the  flat  facet  assumption  the  specular 
term  of  the  BRDF  becomes7 
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where  the  exponent  z  is  given  by 
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where  the  Fresnel  reflectance  is  given  by 

-Reflect  (0  =  Rpec(^)  +  (1  —  Rpec  (^))  (l  —  Usun  ’  u/i)  (4) 

where  Rspec  is  the  specular  reflectance  coefficient.  The  param¬ 
eters  of  the  Phong  model  that  dictate  the  directional  (locally 
horizontal  or  vertical)  distribution  of  the  specular  terms  are 
nu  and  nv .  The  terms  in  Eq.  (2)  are  functions  of  the  reflection 
geometry  which  is  described  in  Figure  2(b).  The  diffuse  term 


(a)  Space  Object  Shape  Model 


(b)  Reflection  Geometry 

Fig.  2:  Reflection  Geometry  and  Space  Object  Shape  Model 


of  the  BRDF  is 

»■(;)  =  a  -  Mi)' 

l_  Un00  'Mobs 

(5) 

where  Rdm(i)  is  the  diffuse  coefficient  for  the  ith  side.  The 
model  discussed  above  assumes  only  single  scattering  and  no 
self  shadowing. 

A.  Flux  Calculation 

The  apparent  magnitude  of  an  SO  is  the  result  of  sunlight 
reflecting  off  of  its  surfaces  along  the  line-of- sight  to  an 
observer.  First,  the  fraction  of  visible  sunlight  that  strikes  an 
object  (and  is  not  absorbed)  is  computed  by 

Run (7)  =  Rsun,vis  (UnW  ’  Usun) 


(6) 


where  Csun,vis  =  1062  W/m2  is  the  power  per  square  meter 
impinging  on  a  given  object  due  to  visible  light  striking  the 
surface.  If  either  the  angle  between  the  surface  normal  and  the 
observer’s  direction  or  the  angle  between  the  surface  normal 
and  Sun  direction  is  greater  than  7r/2  then  there  is  no  light 

reflected  toward  the  observer.  If  this  is  the  case  then  the 

fraction  of  visible  light  is  set  to  Fsun(i)  =  0.  Next,  the  fraction 
of  sunlight  that  strikes  an  object  that  is  reflected  must  be 
computed: 

„  Fsun(i)pm!il(i)A(i)  (u rn(i)  •  u'bs) 

-TohsW  -  ||d/||2  7 

The  reflected  light  of  each  facet  is  now  used  to  compute  the 
total  photon  flux,  which  is  measured  by  an  observer: 

"  N 

F  =  ^  ^obs(^)  +  ^CDD  (8) 

_i=  1 

where  ^cdd  is  the  measurement  noise  associated  with  flux 
measured  by  a  Charge  Coupled  Device  (CCD)  sensor.  The  total 
photon  flux  is  then  used  to  compute  the  apparent  brightness 
magnitude 

raapp  =  -26.7  -  2.51og10  F  (9) 

V2  sun,  vis 

where  —26.7  is  the  apparent  magnitude  of  the  Sun. 

B.  Unscented  Kalman  Filter 

The  attitude  state  errors  are  represented  as  error  Gener¬ 
alized  Rodrigues  Parameters  (GRPs)  resulting  in  a  minimum 
parameter  representation  for  the  attitude  state  error.8  To  within 
first  order,  the  state  error  covariance  of  the  attitude  is  invariant 
whether  the  errors  are  parameterized  using  quaternions  or 
GRPs.9  Therefore,  the  attitude  state  error-covariance  can  be 
directly  decomposed  into  error  GRP  sigma  points  for  use  in  the 
Unscented  Kalman  Filter  (UKF).  The  sigma  points  correspond¬ 
ing  to  the  error  GRPs  are  first  converted  into  error  quaternions 
so  that  the  quaternion  sigma  points  can  be  computed.  The  error 
quaternion,  denoted  by  5q^(i),  associated  with  the  ith  error 
GRP  sigma  point  is  computed  by8 

$Qk  (*)  =  rl  [«  +  5% ,(*)]  XSkP(i)  (10a) 
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where  a  is  a  parameter  from  0  to  1  and  /  is  a  scale  factor, 
which  is  often  set  to  /  =  2  (a  +  1)  so  that  the  attitude  error 
covariance  is  that  of  the  small  roll,  pitch  and  yaw  angle  errors. 
Here  it  is  noted  that  the  subscript  I  and  superscript  B  in  qf 
and  its  estimates  are  omitted  in  this  and  the  following  sections 
for  clarity.  The  ith  quaternion  sigma  point  is  given  by  a  rotation 
of  <5q^  (i)  about  the  a  priori  estimate: 

Ufc  (*)  =  $<lk  (i)  ®  Ufe  (0)  (11) 


The  sigma  points  are  propagated  through  the  system  dynamics: 

X.(i)  =  f  (x(*)>q(*))  (13) 

where 

f  q)  =  4)  —  Jso&j 

v7 
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The  estimated  acceleration  and  torque  due  to  SRP  are  cal¬ 
culated  with  equations  given  in  Ref.  3,  respectively.  After 
propagation,  the  sigma  points  for  the  error  GRP  states  are 
computed  with  the  propagated  attitude  sigma  points.  The 
estimated  mean  sigma  point  quaternion,  q^+1  (0),  is  stored,  and 
error  quaternions  corresponding  to  each  propagated  quaternion 
sigma  point  are  computed  as: 

5qifc+i(*)  =  C+iOO  ®  [qfe+i(0)]_1  (15) 

where  the  notation  for  the  conjugate  quaternion  is  defined  as: 


Using  the  result  of  Eq.  (15),  the  error  GRP  sigma  points  are 
computed  as 


5Pfc+i(*)  =  / 
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After  setting  the  error  GRP  for  the  mean  sigma  point  to  zero, 
the  propagated  sigma  points  are  reconstructed.  The  propagated 
mean  and  covariance  are  calculated  as  a  weighted  sum  of  the 
sigma  points  as 

2  L 

pk+ 1  =  Wi°V  iXk+i(i)  -  x-+1]  [xfc+i(i)  -  X“+1]T 


+  Qk+i 


=  j2wrwxk+i(i) 


where  Qk+i  is  the  discrete-time  process  noise  covariance. 

As  previously  discussed,  measurements  are  available  in  the 
form  of  azimuth,  elevation  and  apparent  brightness  magnitude, 
y  =  [mapp  az  el]T.  Estimated  observations  are  computed  for 
each  sigma  point  using  the  observation  models: 

Vk (i )  =  h  [xk (i) ,  qfc  («)]  (2°) 

The  mean  estimated  output  and  covariance  matrices  are  then 
calculated.  The  quaternion  update  is  performed  by  converting 
the  error  GRP  states  of  to  a  quaternion,  £q^T,  via  Eq.  (10), 
and  performing 

qfc  =  <*qfc  ®  q^  (°)  (2i) 


C.  Multiple  Model  Adaptive  Estimation 


where 


qi  ®q2  =  [’P(qi)  qi]q2 


In  this  section  a  review  of  MMAE  is  shown.  More  details 
(12)  can  be  found  in  Refs.  10  and  11.  Figure  3  shows  the  MMAE 


Fig.  3:  MMAE  Process 


process.  Multiple-model  adaptive  estimation  is  a  recursive 
estimator  that  uses  a  bank  of  filters  that  depend  on  models 
with  different  parameters,  denoted  by  the  vector  p,  which 
is  assumed  to  be  constant  (at  least  throughout  the  interval 
of  adaptation).  Note  the  stationary  assumption  for  the  state 
and/or  output  processes  is  not  necessarily  required  though, 
i.e.  time  varying  state  and  output  matrices  can  be  used.  A 
set  of  distributed  elements  is  generated  from  some  known 
pdf  of  p,  denoted  by  Pr  (p),  to  give  {pW;  £  =  1,  . . . ,  M}. 
The  finite  set  of  parameters  can  be  the  results  of  discretiz¬ 
ing  a  continuous  parameters  space,  selecting  a  set  of  values 
{p^1),  P^2\  •  •  •  5  P^}  dispersed  throughout  the  region  of 
reasonable  parameter  values. 


The  goal  of  the  estimation  process  is  to  determine  the 
conditional  probability  of  the  Ith  element,  p^,  given  all  the 
measurements.  Application  of  Bayes’  rule  yields 


Pr(pW|Yfe) 


Pr(Yfe|pW)  Pr(pW) 

~M 

£Pr(Y*|ptf>)  Pr(p«) 

3  =  1 


(22) 


where  Y/c  denotes  the  sequence  {yo,  yi,  •  •  • ,  y k}-  The  condi¬ 
tional  probability  Pr  (p^  |  Y&)  will  be  the  metric  used  to  select 
the  most  likely  model  and  or  the  most  likely  combination  of 
shape  models.  The  a  posteriori  probabilities  can  be  computed 
through12 


p(yfc|xfc  are  given  as 
Pr(yfe|x“W)  = — - 
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(24) 

where  measurement  residual  for  the  £th  hypothesis  (model)  is 
given  by 

e[e)  =  yk  —  h[x^(pw)]  (25) 

and  corresponding  residual  covariance  matrix  from  the  UKFs 

(26) 


_ pvv  (0 


where  Pkv  ^  is  the  innovation  matrix  for  the  £th  filter. 


Note  that  the  denominator  of  Eq.  (23)  is  just  a  normalizing 
factor.  The  recursion  formula  can  now  be  cast  into  a  set  of 
defined  weights  wk ,  so  that 


=wk- iPr(yfc-i|xfcAi) 
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where  =  Pr(p^|Yfe).  Note  that  only  the  current  time 
measurement  y/c  is  needed  to  update  the  weights.  The  weights 
at  time  to  are  initialized  to  zu^  =  1/M  for  £  =  1,  2,  . . . ,  M. 
The  convergence  properties  of  MMAE  are  shown  in  Ref.  12, 
which  assumes  ergodicity  in  the  proof.  The  ergodicity  as¬ 
sumptions  can  be  relaxed  to  asymptotic  stationarity  and  other 
assumptions  are  even  possible  for  non- stationary  situations.13 


From  Eq.  (24)  and  Eq.  (27)  it  is  seen  that  models  which 
have  lower  residuals  will  have  probability  that  will  increase; 
this  will  favor  models  that  fit  the  observations  better.  Also 
from  Eq.  (24)  it  is  seen  that  models  which  have  small  values 
for  det(S^)  will  have  probability  that  will  grow.  Assuming 
that  all  models  have  same  measurement  noise  covariance 
matrix  R this  will  favor  models  that  have  smaller  variance. 
Therefore  the  MMAE  process  will  tend  to  select  the  maxi¬ 
mum  likelihood  (minimum  variance)  model  from  the  bank  of 
models. 


Pointing 


Projected  Aligned 


Pr(pW|Yfc)  =  Pr(yfc’pW|tfc-l) 

Pr(yfc|Yfc_1) 

=  Pr(yfc|Xfcffl)  Pr(pW|Yfc_i)  (23) 

M 
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The  conditional  probability  of  the  observa¬ 
tions  based  on  each  hypothesis  (likelihood), 


Fig.  4:  Overview  of  Angular  Velocity  Determination  Approach 


D.  Angular  Velocity  Determination 

When  processing  light  curve  observations  it  may  not  be 
valid  to  assume  that  the  SO  is  uncontrolled  and  therefore  we 
must  take  into  account  the  possibility  of  controlled  attitude 
states.  Determining  whether  a  SO  has  active  control  or  not  may 


also  provide  a  feature  state  that  may  be  used  for  classification. 
For  example,  a  determination  of  whether  a  SO  is  passive  or 
active  can  be  made  based  on  whether  light  curve  observations 
indicated  that  the  SO  has  active  attitude  control. 


In  this  work,  the  attitude  control  is  simulated  by  assuming 
control  profiles,  for  example  Sun  pointing,  Nadir  pointing,  and 
spin  stabilized.  Then  for  each  control  profile  a  desired  angular 
velocity  is  determined  which  will  allow  the  SO  to  track  the 
relevant  directions.  The  angular  velocity  profiles  are  used  to 
calculate  the  torque  required  to  track  this  profile.  This  section 
discusses  the  attitude  control  approach  used  and  four  method 
for  calculating  the  desired  angular  velocity  profile  which  is 
shown  in  Figure  4.  The  attitude  control  is  designed  to  minimize 
the  following  error: 

e  =  a?  -  ud  (28) 

Differentiating  this  equation  with  respect  to  time  yields 

e  =  dj  -  ud  (29) 

It  is  desirable  for  the  error  dynamics  to  decay  exponential  over 
time,  i.e.  e  oc  e~kpt ,  and  therefore  the  error  rate  equation  is 
desired  to  have  the  following  form: 

e  =  —  kpe  (30) 

Then  using  Euler’s  equation  and  assuming  disturbance  torques 
are  negotiable,  Eq.  (29)  can  be  written  as 

e  =  Jsq1  (t  -  [ux]  Jso^)  -  Ud  (31) 

where  r  is  the  torque  provide  by  the  attitude  actuator.  Then 
for  an  exponentially  decaying  tracking  error  the  desired  torque 
expression  becomes 

r  =  ~Jso  (T  -  [u?x]  Jsov)  +  Vd  -  kpe  (32) 

This  expression  is  used  to  calculate  the  torque  required  to 
maintain  the  desired  pointing  profile. 


In  this  section  angular  velocity  determination  approaches 
are  discussed.  Consider  the  following  unit- vector  measurement 
model  at  time  tk: 

bjk  =  Ak  rj  +  vifc  (33) 

where  b7fc  is  the  jlh  pointing  vector  in  the  inertia  frame  and 
is  Yj  the  same  pointing  vector  in  the  body  frame.  The  attitude 
matrix  mapping  from  inertia  to  the  body  frame  is  denoted  by 
Ak .  Our  goal  is  to  determine  the  rate  of  change  of  this  attitude 
matrix  or  the  angular  velocity.  Taking  the  difference  between 
successive  measurements  of  Eq.  (33)  gives 

bjk+1  -  b Jh  =  [Ak+ 1  -  Ak]  r ,•  +  vifc+1  -  wjk  (34) 

We  assume  that  the  body  angular  velocity  u  is  constant 
between  tk  and  tj. c+i,  and  ignore  terms  higher  than  first  order 
in  u>  At.  With  these  assumptions  the  following  first-order 
approximation  can  be  used:14 


Ak+i 


hx3 


At  [ukx] 


Ak 


(35) 


In  this  case  is  the  average  velocity,  but  this  becomes  less 
of  a  problem  as  the  sampling  interval  decreases.  Substituting 
Eq.  (35)  into  Eq.  (34)  gives 

b jk+1  ~  bjk  =  -At  [wfe x]  Ak  rj  +  vjk+1  -  vjk  (36) 


Our  goal  is  to  determine  an  angular  velocity  estimate  indepen¬ 
dent  of  attitude  and  the  reference  vectors.  This  is  accomplished 
by  solving  Eq.  (33)  in  terms  of  Ak  and  substituting  the 
resultant  into  Eq.  (36),  which  yields 

^  [b jk+1  ~  bjk }  =  [bjk  x  ]  cjfe  +  wjk  (37) 


where  w jk  is  the  new  effective  measurement  noise  vector  given 
by  , 

wj,  s  -  vj.l  (38) 

Note  that  At  will  have  finite  values,  since  discrete-time  mea¬ 
surements  are  assumed.  Equation  (37)  can  now  be  cast  into  a 
linear  least- squares  form  for  all  measurement  vectors,  which 
leads  to 


nk 

l^[bj*x]T-R7h1[bJ-fcx] 

-i=i 


3»XfRju1(*>3H+l 


3=1 


biJ 


(39) 


where  uk  is  the  estimate  of  For  small  At  the  propagated 
true  value  of  b j  can  be  given  using  Eq.  (35): 

bjfc+1  ~  {/3X3  -  At  [wfcx]}  bifc  (40) 

Substituting  Eq.  (40)  into  Eq.  (35),  left  multiplying  by  [bJfc  x]T 
and  right  multiplying  by  [b jk  x  ]  gives 

[bjk  x]rJR“1[bife  x]  =  a~2[bjk  x]T[bjfe  x]  (41) 

where  =  2cr|/Af2.  Also,  since  bJfc+ibJfc  «  1,  it  is  easy  to 
show  that 

[bjkx]TRjk{bjk+ 1  ~bjk)  «  vJ2[bjkx]Tbjk+1.  Therefore, 
Eq.  (39)  is  well  approximated  by 
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jk  +  1 


(42) 

where  the  measurements  have  again  been  substituted  in  place 
of  their  true  values. 


III.  Simulation  Results 

Four  simulation  scenarios  are  presented  to  show  the  perfor¬ 
mance  of  the  MMAE  based  classification  approach  to  classify 
an  SO  from  magnitude  and  angles  observations.  In  each 
scenario  a  different  object  is  selected  that  falls  into  a  different 
class.  The  objects  selected  are  spin  stabilized  bus,  uncontrolled 
bus,  nadir  pointing  bus,  and  uncontrolled  rocket  body.  All 
scenarios,  an  SO  is  in  near  geosynchronous  orbit  with  orbital 
elements  given  by  a  =  42,364.17  km,  e  =  2.429  x  10-4, 
i  =  30  deg,  cu  =  Q  =  0.0  deg  and  M0  =  91.065  deg.  The 
simulation  epoch  is  15 -March-20 10  at  04:00:00  GST.  The 
initial  quaternion  and  angular  rate  of  the  SO  are  given  by 
qf  =  [0.7041  0.0199  0.0896  0.7041]r  and  u3§/f  = 

[206.26  103.13  540.41]T  deg/hr. 

Brightness  magnitude  and  angle  observations  are  simulated 
using  a  ground  station  located  at  20.71°  North,  156.26° 
West  longitude  and  3,058.6  m  altitude.  Measurements  con¬ 
structed  using  instantaneous  geometry  are  corrupted  by  zero- 


mean  Gaussian  white  noise  with  standard  deviations  of  1 
arc-seconds  on  the  azimuth  observation,  1  arc-seconds  on 
the  elevation  observation  and  0.1  for  the  brightness  magni¬ 
tude.15  Observations  are  available  every  5  seconds  for  one 
hour. The  initial  states  for  each  filter  are  given  by  qf  (to)  = 
[0.7500  0.0712  0.0947  0.6508]T  (a  10  degree  attitude 

error),  u>f/7(t0)  =  [220.26  117.13  554.41]T,  a(t0)  = 

42,364.148255  km,  e(t0)  =  2.4290  x  1(T4,  i(t0)  =  30.0083 
deg,  u)(to)  =  —1.172  deg,  Cl(to)  =  0.0  deg  and  Mo  (to)  = 
92.137  deg.  Initial  3cr  values  are  taken  to  be  20  deg  for 
the  attitude  states,  72  (deg/hr)  on  the  angular  rates,  300  km 
on  position  and  3  (km/s)  on  velocity.  The  process  noise  for 
the  UKFs  are  taken  as  QKk  }  =  0  for  this  proof  of  concept 
simulation. 


A.  Classification  Results 

Figure  5  shows  the  classification  results  for  a  spin  stabi¬ 
lized  bus.  In  this  case  the  bus  models  are  considered  to  be 
regular  cuboids  with  aspect  ratio  larger  than  0.1.  As  discussed 
in  classification  section  there  are  a  number  of  classes  the 
classification  approach  determines,  whether  the  SO  belongs 
to  these  classes.  The  first  class  is  whether  the  SO  is  intact  or  a 
fragment.  In  this  case  the  true  model  is  an  intact  bus  and  from 
Figure  5  we  can  see  that  this  determination  is  made  relatively 
quickly.  The  second  determination  is  whether  SO  is  active  or 
passive,  which  is  also  shown  in  Figure  5.  The  active  or  passive 
decision  is  made  using  the  probability  of  all  active  and  passive 
models  in  the  bank.  We  can  see  from  Figure  5  that  the  active 
model  gain  achieves  large  probabilities  after  0.2  hours. 

Additional  examples  are  shown  for  an  uncontrolled  bus 
(Figure  6),  nadir  pointing  bus  (Figure  7),  and  uncontrolled 
rocket  body  (Figure  8).  The  approach  shows  good  performance 
for  these  examples  for  determining  the  correct  class.  From 
the  figures  it  can  be  seen  that  some  objects  take  longer  to 
classify.  This  is  due  to  the  fact  that  for  some  spin  states  the 
light  curves  are  similar,  but  for  uncontrolled  spin  states  the 
light  curve  differs  significantly.  This  can  be  seen  from  nadir 
pointing  bus  (Figure  7),  which  takes  the  longest  to  converge 
to  its  classification,  and  from  uncontrolled  rocket  body  (Figure 
8),  which  converges  the  fastest  to  its  classification. 


(a)  Class  1:  SO  Size  Fea¬ 
tures 


(b)  Class  2:  Attitude  Con¬ 
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(c)  Class  3:  SO  Type 


(d)  Class  4:  Spin  Control 
State 


Fig.  5:  Spin  stabilized  Bus  Example 
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IV.  Conclusion 

In  this  paper,  an  Multiple  Model  Adaptive  Estimation 
(MMAE)  scheme  for  space  object  classification  using  light 
curve  and  angles  data  was  presented,  which  can  be  used 
to  identify  the  most  probable  class  of  the  SO  along  with 
its  associated  rotational  and  translational  states.  Using  an 
Unscented  filter  to  reduce  brightness  magnitude  and  angle  data, 
the  MMAE  is  able  to  determine  the  probability  of  each  model. 
An  approach  is  presented  that  uses  the  probability  from  a 
MMAE  process  to  determine  the  probability  that  a  given  Space 
Object  (SO)  falls  in  a  given  class.  The  classification  approach 
determines  whether  the  SO  is  intact  or  fragment,  its  control 
states,  the  type  of  control  state,  and  whether  it  is  rocket  body, 
payload,  or  debris.  Simulation  results  showed  a  number  of 
examples  and  good  results  for  classification  are  given. 


(c)  Class  3:  SO  Type 


(d)  Class  4:  Spin  Control 
State 


Fig.  6:  Uncontrolled  Bus  Example 
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Fig.  7:  Nadir  Pointing  Bus  Example 
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(c)  Class  3:  SO  Type 


(d)  Class  4:  Spin  Control 
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Fig.  8:  Rocket  Body  Example 


